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Abstract 

Cortical neurons include many sub-cellular processes, operating at multiple timescales, 
which may affect their response to spiking stimulation through non-linear and stochastic 
interactions with ion channels and ionic concentrations. Since new processes are con- 
stantly being discovered, biophysical neuron models increasingly become "too complex to 
be useful" yet "too simple to be realistic". A fundamental open question in theoretical neu- 
roscience pertains to how this deadlock may be resolved. In order to tackle this problem, 
we first need to define the notion of an "excitable neuron model". Then we analytically 
derive the spiking Input-Output (I/O) relations of such neuronal models, relating input 
spike trains to output spikes based on known biophysical properties. We show that the 
mean firing rate I/O is a memoryless non-linear function, while the relation between the 
input spike-times to the output spike-occurrences is given by a linear stochastic dynamical 
system, with parameters tuned by the mean stimulation rate. Thus we obtain closed-form 
expressions for the mean firing rates, all second order statistics (input-state-output cor- 
relations and spectra) and construct optimal linear estimators for the neuronal response 
and internal state. These results are guaranteed to hold, given a few assumptions (mainly, 
a timescale separation assumption), for any stochastic biophysical neuron model (with an 
arbitrary number of slow kinetic processes) under general sparse stimulation. We numeri- 
cally test the resulting expressions for various models, and show that they hold well, even 
in cases our assumptions fail to hold. Our results suggests that the common simplifying 
approach that ignores much of the complexity of the neuron might actually be some- 
times unnecessary and even deleterious. Specifically, the stochasticity of ion channels and 
the temporal sparseness ("spikiness") of the inputs is exactly what rendered our analysis 
tractable, allowing us to incorporate multiple slow kinetic processes. 



Introduction 

When aiming to describe a complex system mathematically, it is important to use models 
that are both simple enough to understand yet complex enough to capture reality at a 



required level of detail. This tradeoff seems to be ubiquitous in quantitative biology in 
general [I], and in neuroscience in particular |2j. Since the seminal work of Hodgkin 
and Huxley [3] on the squid's giant axon, biophysical models of neurons have become 
an established framework through which many experimental results can be contrasted, 
reconciled, interpreted and related to functionality (e.g., (1] and references therein). Such 
models, termed "Conductance Based Models" (CBMs), are commonly used to investigate 
intra-cellular dynamics (e.g., dendritic integration |5|, axonal propagation (6] and ionic 
currents [7]) and their relations to neuronal firing patterns |8}[9]. 



However, CBMs are less popular when modeling networks of neurons 10-14] or even 
when attempting to accurately predict the response of a single neuron to stimuli (2]. 
The main reason for this is that CBMs are complex non-linear models, containing many 



variables and unknown parameters (sometimes ranging in the hundreds 15 16)), not all 
of which can be identified [17]. Moreover, the response of CBMs can be highly sensi- 
tive to values of these parameters - precisely because CBMs model excitable cells. These 
properties usually render CBMs notoriously difficult to tune |2], highly susceptible to over- 
fitting [2], computationally expensive |18|, mathematically intractable [2] and generally 
"hard to understand" - in the sense that given a specific stimulation pattern and parameter 
values, it is hard to quantitatively, and even qualitatively, describe the neuronal response 
of the CBM. Therefore, modelers often prefer to use simplified, phenomenological or sta- 



tistical versions of CBMs 19 -261, thereby reducing the biophysical interpretability of the 



model. Moreover, it is not clear when such approaches neglect essential biophysical con- 
straints or overly simplify the neural response so that it looses some important functional 
attributes. 

To make matters worse, it now seems that even very detailed CBMs should be strictly 
considered as highly simplified models. New sub-cellular kinetic processes are being dis- 



covered at an explosive rate 27 29 . This variety is particularly large for slow kinetic 
processes 30 which may affect the neuronal response, as indicated in a recent experi- 
ment (3l]. Additionally, most CBMs tend to ignore the stochastic nature of ion chan- 
nels |32[|33| , which was shown to affect the neuronal response and functionality ( |3~4~|[3~5"] 
and references therein). It would be natural to believe that incorporating into the CBM 
framework many additional biophysical processes as well as stochasticity would lead to 
intractable and unwieldy equations. Surprisingly, as we show here, the opposite is the case 
when temporally sparse inputs are used (so that they are "spiky" - short and well-spaced). 
In that case, we show that stochastic CBMs can become significantly more tractable, in- 
dependently of the number and complexity of the slow kinetic processes added. Before we 
go into details about our results and methodology, we explain why a sparse stimulation 
regime is physiologically reasonable. 

It is usually assumed that the synaptic inputs to cortical neurons are rather homo- 
geneous and uncorrelated, so that the neuronal firing is dominated by the summation of 
many (~ 10 4 ) such similar inputs ("the diffusion approximation" - see (36J and references 
therein). However, recent evidence suggests that this is not generally true. First, the mean 
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firing rates of cortical neurons are rather low - less than 1Hz ( |37|; also 38 and references 
therein). Furthermore, incoming synaptic inputs are usually integrated sequentially and 
non-linearly at several levels of dendritic compartments (39|-|4T], where adjacent synaptic 



spines are frequently synchronized |42|. Moreover, recent studies have shown that the 



distribution of synaptic efficacies in the cortex is log-normal 43 , 44] - so a few synapses 
are very strong, while most are very weak. This indicates that the neuronal firing patterns 
might in fact be dominated by a small number of very strong synapses while the sum of 
the weak synapses sets the voltage baseline (45). Such a possibility is supported by the 
fact that stimulation of a single cortical cell can affect the network activity with a single 



spike |46|, elicit whisker movements in rats |47| and even modify the global brain state 48 



Taken together, these observations suggests that the above-threshold stimulation reach- 
ing the soma may be temporally sparse, as can be observed directly in some cases (36] . 
Note that in other cases it might be difficult to observe this due to sub-threshold voltage 
fluctuations caused by weak synaptic inputs or intrinsic neuronal dynamics. 

Now suppose, as reasoned above, that incoming inputs to a stochastic CBM are indeed 
temporally sparse (Fig. [T]A.) - so Action Potential (AP) duration is "short" in comparison 
to the inter-stimulus interval (Fig. [Tj3) . In that case, we can construct a discrete time 
stochastic map describing the dynamics of the neuronal state and response, sampled at 
each stimulation time (Fig. [Tp, top). This map can be accurately and transparently de- 
scribed using parameters directly derived from the original CBM, if we additionally assume 
that the slow kinetics modulating neuronal excitability are "slow" in comparison to the 
inter-stimulus intervals. The resulting "excitability map" (Eqs. 20 21), has a significantly 
smaller number of parameters than the original CBM, while maintaining a well-defined 
biophysical interpretation. Using this map, simulation time can be dramatically reduced, 
using a stimulation-event based simulation. 

Furthermore, under stationary input, if we further assume that AP generation dy- 
namics is sufficiently "noisy", the excitability of a cell will fluctuate near a fixed point. 
Linearizing the excitability map around this "excitability fixed point" we derive explicit 
mathematical expressions for the neuronal input-output relation. Specifically, we find 
that the mean firing rate can be written as a function of the mean stimulation rate, while 
the fluctuations can be described using a stochastic dynamical linear system, with param- 
eters tuned by the mean stimulation rate. This linear system yields a highly interpretable 
biophysical description of the quantitative dynamics (Fig. [Tj bottom), which contains a 
feedback loop, similarly to many other biological systems (e.g., |49|). This linear descrip- 
tion allows us to use the well established engineering mathematical toolbox to analyze 
the system, write explicit expressions for the second-order statistics of the neuron (all 
input-state-output correlations and spectra) and develop linear optimal estimators for 
the neuronal state and output. Lastly, we show numerically that our approximations hold 
well in many cases, even when some of our assumptions (e.g., timescale separation and 
fixed point assumptions) break down. 

In |35| and more recent unpublished work we demonstrate the applicability of our 
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approach to the experimental data of 31 , where synaptically isolated individual neurons 
(from rat cortical culture) were stimulated with physiologically sparse inputs for extended 
durations of days. Interestingly, the neurons responded in a complex and irregular man- 
ner ("1// spectrum) over the entire range of experimental timescales - from seconds to 
days. Due to the large number of processes involved at such timescales, many of them un- 
known or lacking known parameters, this phenomenon cannot be modeled using a purely 
simulation-based approach. However, using the analytic tools developed here we can cir- 
cumvent this problem, exactly characterize the classes of CBMs that can generate the 
observed behavior (given our assumptions), and perform parameter estimation from the 
observed neuronal response. 

Similarly, our method could be immediately applicable in other cases where a direct 



pulse-like stimulation of a neuron is used, as an experimental probe (e.g., 50,51]) or 
in implants (e.g., in the cochlea; see |52| and references therein). More generally, we 
hope that the generic methods described here (or some variant of them) could be used 
to significantly improve simulation, analysis and design of CBM-based networks, while 
allowing CBMs' complexity to increase according to experimental findings. For example, 
similar methods may help provide a definitive characterization of the effects of variability 
|53~] and co- variability [54] of neuronal parameters on its response in some cases, since we 
have closed-form expressions relating neuronal response to biophysical parameters. 

The remainder of the paper is organized as follows. We begin in the 'Models' section 
by presenting the general mathematical setup of an abstract neuron under temporally 
sparse stimulation ; and define formally the general class of CBMs. Then, in the 'Results' 
section we show how to reduce such systems to simple discrete time stochastic maps; 
derive simple, yet precise, input/output relations for the neuron at the level of spike 
trains; derive input-output-state second order statistics; explain how to perform linear- 
optimal estimation of internal neuronal variables and predict firing patterns; and finally, 
we demonstrate our results numerically for various biophysical models. Recall, however, 
that our main purpose here was to derive very general results - which are valid for a broad 
class of biophysically realistic conductance based models. These results are summarized 
and set in context in the 'Discussion' section. In Information SI we delve deeper into 
the mathematical derivations and expands on some of the numerical issues involved. In 
Information S2 we give the details of the various model parameters we used in simulations. 
Finally, in Information S3 we list all of our mathematical notations, for convenience. 



Models 

General formalism 

Problem definition Isolated cortical neurons are excitable elements that can generate, if 
stimulated, an Action Potential (AP) response - a "spike" in the cell's membrane voltage. 
Formally, the "AP event" is defined by a some specific condition. For example: "an AP 
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has occurred if V, the voltage at the axon hillock, has crossed upwards some threshold 
(say, — 10mV)". As seen experimentally [3l] , APs are not generated spontaneously by 
cortical neurons if no stimulation is given (either synaptically or directly) . Therefore, we 
assume that this is the case (i.e., we do not consider oscillatory neurons). Additionally, 
in this work we consider short temporally sparse stimulation events. These stimulations 
can represent either stimulation by direct current pulses (as in |3l]), single pre-synaptic 
APs, or simultaneously occurring groups of pre-synaptic APs. We are interested in finding 
the "spiking" input-output relation for these neurons under sequences of such stimulation 
events. 

Formally, suppose that at times {t m }^ =0 a neuron receives a sequence train of identical 
stimuli. We denote by {1^}^ =0 the occurrence events of APs at times {t m }^ =0 , i.e., after 
each stimulation time t m (Fig. [T]A.) 



Defining T m = t m+ i — t mj the inter-stimulus interval, and r^p as the upper timescale of 
an AP event (Fig. [!£) we assume that 

1. The neuron does not spontaneously fire APs. 

2. The spike times {t m }^ =0 are temporally sparse, i.e. tap <C T m for every m ("no 



Objective: Mathematically characterize the relation between {Y m } and {T m } under the 
most general conditions. 




1 , if an AP occurs 
, otherwise 



(1) 



collisions"). 
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Figure 1. Mathematical analysis - schematic summary: A Main objective - find 
the input/output relation between inter-stimulus intervals (T m ) and Action Potential 
(AP) occurrence events (Y m ) - for any biophysical neuron model. B An AP is said to 
have occurred if the voltage V crossed a threshold Vth following the stimulus. We 
examine a "sparse" stimulation regime in which T m ^> tap, the timescale of AP 
generation (so stimulation pulses do not "collide"). C Model reduction scheme. Each 
number references the relevant (sufficient) assumptions in the text. The end result is the 
conversion of a complex biophysical neuron model to a simple linear model (Eqs. 27 28), 
with parameters (F, d, a...) directly linked to the biophysical parameters. 
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Probabilistic approach Since the specific details of internal neuronal dynamics are 
not fully understood, we will attempt to answer this question with as few assumptions 
as possible, in order to maintain a high level of generality. We begin with a rather 
abstract neuronal model, in which x (t) is a vector representing the internal neuronal 
"state" at time t. For example, in the context of the Hodgkin-Huxley model |3| x (t) 
represents the voltage, the sodium channel variables and the potassium channel variable. 
We define x m = x (t m ) and the history set H m = {{^kYYo , {Tfc}£l > {**}fclo} ( note tnat 
Hm C Tim+i)- Assuming that x m is a "state vector" implies the Markov property, namely 
it is a sufficient statistic on the history to determine the probability of generating an AP 
at each stimulation, 

PAP ( x m) = P (Ym = l| x m) = P (Ym = 1 | x m.) T^m-l) i (2) 

and, together with Y m and T m , its own dynamics 

P ( x m+l | x m; T m , Y-m) = P ( x m+l \H-m) i (3) 

which implies the following causality relations 

Tn Tl 

xo 4 xi 4 x 2 ••■ x m 

I / I / I I • ( 4 ) 

Yo Yi Y 2 ■ ■ ■ Y m 



This causality structure is reminiscent of the well known Hidden Markov Model |55|, 
except that in the present case the output Y m , affects the transition probability, and we 
have input T m . 

Theoretically, if we knew the distributions in [2] and |3j as well as the initial condition 
P (xo), we could integrate and the find an exact probabilistic I/O relation P ({Yk}™ =0 \ {^fc}fcL ) 
However, it may be hard to find an expression for P (x m+1 |x m , T m , Y m ) in general bio- 
physical models, and also the integration might be intractable. We therefore propose an 
alternative approach in the next section. 

Stochastic map approach We construct a dynamic "state-space" system with T m , the 
inter-stimulus interval lengths, serving as inputs, x m representing the neuronal state, and 
Y m the output. The model will be augmented with two noise sources. The first noise 
source is the "innovation part" |56| of the AP generation process, 

Y m (Y^Jxjjj, T~L m ^\) Y m Pap ( x m) > (5) 

where (X\Y) is the conditional mean of X given Y. The second noise source is the 
innovation part for the state increments Ax m = x m+ i — x m 

n?n ^ x m (^ x m | x m; -^mj ^m) • (6) 
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Interestingly, we can prove (see Information SI, section A) that e m and n m are un- 
corrected white noise processes - i.e. (e m ) = 0, (n m ) = (e n n m ) = and (n m n^) = 
(n m n^) 5 mn , (e m e n ) = (e 2 m ) 5 mn where 5 nm = 1 only if n = m and otherwise. Now, 
defining h + (x m ,T m ) = (Ax m |x m , T m , Y m = 1) and h (x m ,T m ) = (Ax m |x m , T m , Y m = 0) 
and using Eqs. [5] and [6] it is straightforward to obtain the following non-linear, stochastic, 
dynamical system 

^m+l -|- l^^ll (x m , T m ) -|- (1 Y m ) h (x m , T m ) -\- Il m j (7) 

= Pap (x m ) + e m . (8) 

Formally, such a system could be used to simulate the neuronal response. Note that in 
such a simulation, pap (•) and h (•, •) are deterministic functions while n m and e m are 
random processes conditioned on (x. m ,Y m ,T m ) and x m , respectively. In general, in order 
to perform this simulation, we need to know pap (•)> h ± (■, ■) and the full distribution of 
n m given (x m , Y m ,T m ). Although both the probabilistic and stochastic map approaches 
are formally equivalent, the latter approach can be more useful in some cases, as we shall 
later see. For concreteness we examine the canonical class of biophysical neuron models 
- Conductance-Based models (CBM). 



Stochastic conductance-based models 

Due to our objective, we are currently not interested in modeling how stimulations from 
different neurons are integrated in the dendrites or how the propagation in the axon affect 
the latency of the AP. Therefore, in order simplify the following definitions and analysis, 
we avoid unnecessary generality, and consider only models with a single (isopotential) 
compartment. In a general stochastic CBM for such a point neuron the voltage dynamics 
is determined by ion channels - protein pores that open and close stochastically in a voltage 
dependent way |33|. To formalize this, the dynamics of each channel is described as an 
irreducible continuous time Markov process (57| , with voltage dependent rates. Given the 
voltage, all channels are assumed to be independent. Such a mathematical framework 
is also able, in many cases, to include other cellular processes affecting excitability (e.g., 
changes in ionic concentrations). However, for descriptive simplicity, we will refer to these 
Markov processes as "ion channels" only. We index by c the different types channels, 
c = 1, C. For each channel type c there exist iV( c ) channels, where each channel of type 
c possesses K^ c ' internal states. Also, the kinetic transition rate from state j to state i in 
a c-type channel is given by A\j (V). And so, for each channel of type c that resides in 
state i, the probability that the channel will be in state j after time dt is given by 

V-Z^^(v)dt , if i = j • u 

We define to be the fraction of c-type channels in state k, x.^ to be a column 
vector composed of c = 1,2, ... ,C, and x to be the vector formed by concatenating 
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the x( c ) vectors. The dynamics of the voltage V is then determined by V, x and / (t), the 
external current stimulation, 

v = f(y,x,i(t)) , (io) 

for an appropriate function /. Note that the neuronal state is now (V, x), and not x as 
was before. For an excitable point neuron model, rather than some general dynamical 

system, both / and |y4-^| must have adhere to additional constraints. 
The excitability constraint 

What is an "excitable neuron"? The answer to this question must be formalized math- 
ematically before any generic analytical results can be derived regarding such elements. 
The main characteristic of the classic HH model |3i| is the almost binary "all-or-none" 
"short-duration" response of its variables to pulse stimulation. It is desirable that an "ex- 
citable element" will also have this "all-or-none" "short" response. However, if, as in later 
models, some of the parameters of the model are changed into slow variables (e.g., slow 
sodium inactivation [58]), then the response of these variables is neither "all-or-none" nor 
"short". Can we extend the definition of "excitable" to such cases? 

To be able to do this, we need to assume that we can express x as a combination of a 
"rapid" component r, which generates the shape of the AP (together with V), and a "slow" 
component s, which modulates the system over larger timescales. Note that both r and s 
are vectors, comprising multiple components. For example, if only rapid or slow channel 
types exist, then x is simply the concatenation of r and s. For simplicity of notation, we 
shall henceforth assume that this is the case. However, such a decomposition is readily 
extendable to cases where some channel types possess both slow and fast kinetics (as 
we later demonstrate; see also (59)). And so, if tap and r s are the respective kinetic 
timescales of {V, r} and s, then tap < r s . Using this decomposition, suppose we "freeze" 
the dynamics of s (so that effectively r s =oo) and allow only V and r to evolve in time. We 
can say the original model describes an "excitable point neuron", if the following conditions 
hold in the "semi-frozen" model: 

1. If / (i) = 0, then for all initial conditions, V and r rapidly (within timescale tap) 
relax to a constant and unique steady state ("rest"). 

2. Assume that V and r are near rest, and a short stimulation pulse is given with 
duration t stim < tap and amplitude / s tim- For certain initial conditions and val- 
ues of Jstinx, we get either a stereotypical "strong" response ("AP response") or a 
stereotypical "weak" response in V ("no AP response"). Only for a very small set of 
initial conditions and values of J s ti m , do we get an "intermediate" response ("weak 
AP response"). By "stereotypical" we mean that the shape of response does not 
change much between trials or for different initial conditions in {V, r} (however, it 
can change with s). 
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Note this definition for an excitable element can be easily extended to be used even if the 
model is not strictly a stochastic CBM. 

Stochastic differential equations form The original state vector x (from section 
'General formalism') was split in the previous section into three sets of variables: V 
(voltage), r (rapid kinetic variables) and s (slow kinetic variables). The variables V and 
r determine AP generation, and so are characterized by the upper timescale of tap- The 
"excitability" variable s modulates the AP generation process and operates at a lower 
timescale of r s . Using the method similar to that derived in (60) (see Information SI, 
section B) we can write an alternative description for this model using the following 
stochastic differential equations 



V = f(V,T,S,I(t)) (11) 

r = A'(V)r-b'(V)+B'(V,r)g (12) 
s = A(V)s-b(V) + B(V,s)£ (13) 

where / is some function, £and£' are vectors of independent white noise processes with 
zero mean and unit variance, A, A', B, B', b and b are given explicitly as functions of 
the microscopic voltage-dependent kinetic rates Af, , and B, B' also depend on s, r and 
the numbers of channels N^ c \ Using the above equations the dynamics of a stochastic 
CBM can be simulated efficiently (for details see Information SI, section C). Note both 
A and A' are stable matrices and D =BB T and D'=B'B /T are both positive definite. 
Therefore, if V is held constant, the (ensemble) means (s) and (r) te nd to S qq = A x b 
and Too = (A / ) _1 b / , respectively. Note that the linear nature of Eq. 
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13 



stems from 

the fact that ion channels are assumed to be independent given the voltage. Also, the 
fact that these equations are not directly coupled is not essential and is made only for 
mathematical convenience (our results hold for the coupled case, as can be seen in Fig. 



Results 

The excitability map 



Given a sequence {£ m }^ =0 of stimulations (as defined in section 'Problem definition'), we 
would like to find explicitly how V m = V (t m ), r m = r(t m ) and s m = s (t m ) in a CBM 
(section 'Stochastic conductance based models') evolve in time, as in Eqs. [7]|8] (with 
x m = (V m , r^, s^) T ). Theoretically, this could be done by a direct integration of the 



CBM Eqs.(ll 13J) between stimulations. However, such integration does not generally 
admit an interpretable closed form expression. To overcome this, we will need to assume 
that the slow kinetics are significantly slower than the inter-stimulus intervals 
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Assumption 1: T m <C r s for every m. 

This assumption, together with our assumption on the stimulus (tap <C T m for every m, 
with tap possibly extended so to include the timescale of r) and the excitable nature of the 



CBM, renders the dynamics (11 13) between stimulations relatively easy to understand. 
Specifically, between two consecutive stimulations, the slow variable s (£), starting from its 
initial condition at the time of the previous stimulation, changes according to the voltage 
trace V (t) (which enters through the parameters in Eq. 13). At the same time, the 
fast variables (V (t) , r (t)) follow stereotypically either the "AP response" (Y m = 1) or the 
"no-AP response" (Y m = 0), then equilibrate rapidly to some quasi-stationary distribution 
q(V,r\s). Summarizing this mathematically, we obtain the following approximations 

P (^mJSm) ~ j P (Y m \V m , F m , S m ) (J (Vrni T m \s m ) dV m dY m , (14) 

Using these equations together with Eqs. [21 [3j we obtain 



Pap ( s m) — P (Y m — l|s m ) — P (Y m — l|s m , Hm-i) > (16) 

P (Sm+l|Sm) ^mj I'm) ^ (Sm+iPmi ^mi -^mi ^m-l) > (I'O 

where, slightly abusing notation, we redefined pap (■)• Therefore, due to our assumptions, 
the "excitability" vector s m can now replace the full state vector x m = (V m , r 7 ^, s^) as the 
sufficient statistic that retains all relevant the information about the history of previous 
stimuli. Picturesquely, this gives us a Markov process with the causality structure 

So — > Si —7- S2 • ■ ■ S m 

I / I / I ; ■ (is) 

Yq Y\ i2 • • ■ Y m 
Since the function j?ap ( s ) is not affected by the kinetics of s, it can be found by numerical 



simulation of a single AP using only Eqs. 11 12, when s is held constant (see Information 
SI, section D). 

Next, we aim to find a recursive stochastic map for s m , similarly to Eq. [7j Integrating 



Eq. 13, we obtain 



m+1 



s m+ i = s m + / (A(V(t))s(t)-b(V(t)) + B(V(t),s(t))£(t))dt (19) 

where A, b, B and D = BB T depend on the voltage through the kinetic rates. To calculate 
this integral, let a (V) be some microscopic voltage-dependent kinetic rate (some (V), 
as defined in section 'Stochastic conductance based models'), and average it between two 
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consecutive stimulation pulses given Y mi s m and T m . We obtain 



a{V (t))dt\s m ,Y m ,T m 



tap 



1 

tap 



-TAP 



(a+ (s m ) Y m 



a +(T m -r AP ; 
;s m ) (1 - F m )) r AP + (T m - 



T m — tap 
tap) «o (s m ) 



a{V{t))dt 



tm+7 AP 



where in the last line we defined ao,a + and a_ to represent the average rate at rest, 
during an AP response and during a no-AP response, receptively. Similarly, we augment 
A, D, B and b with {0, ±} subscripts to denote the same quantities, when all the kinetic 
rates are substituted with their averaged {0, ±} equivalents. Since A, D and b are linear 
combinations of the kinetic rates (see Information SI, section B), this gives 

tm + 1 

A (V (*)) dt\Y m , s m , T m = (A+ (s m ) Y rn + A_ (s m ) (1 - Y m )) r AP +(T m - r AP ) A (s m 



and similarly for D and b. Using this notation, we can perform the integration in Eq. 19 



we obtain a similar dynamical system as in Eqs. [7] and [8] (in which x m is replaced with 
s m ) and obtain the "excitability map" for the neuronal dynamics 

+ tap (Y m (A + (s b+ (s m )) + (1 - Y m ) (A_ (s m ) s m - b_ (s 



'm+l 



Y, 

with 
(n m n 



+ (T m - r AP ) (A 
= Pap (s m ) + e m . 



(s m )) + n r 



Y m ) w tap (l" m D + (s m ) + (1 - Y m ) D_ (s m )) + (T„ 



Second order corrections in the mean for Eq. 



Eq. 



22 



20 



are of order O ^(T m r t 



tap) Do (s m 
-i\ 2 



(21) 

• (22) 
while 



has corrections of size O ((TrnT- 1 ) 2 /JV") where N = min c A^^ (N^ is the channel 

number of the c-type channel, defined in section 'Stochastic conductance based models'). 
Henceforth we shall neglect these corrections, since r m r s _1 <C 1, by assumption 1. 

Additionally, the distribution of n m given s m , T m , Y m can be generally computed using 



the approach described in |60|. For example, it can be well approximated to have a 
normal distribution if channel numbers are sufficiently high and channel kinetics are not 
too slow 60 . In that case only knowledge of the variance (Eq. 22 ) is sufficient to generate 

Urn • 



And so, using Eqs. 16, 20 and the full distribution of n m , we can now simulate the 
neuronal response using a reduced model, more efficiently and concisely (fewer parameters) 
than the full CBM (Eqs. [Tlp3|, since every time step is a stimulation event. Note that 
the reduced model parameters, having been deduced from the CBM itself, still retain a 
biophysical interpretation. Next, we wish to further simplify this reduced model, in the 
case of stationary input. 
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Input-output relations for a stationary input 

We now examine the case where {T m } is a Wide Sense Stationary (WSS) process, with 
mean T*, so that the assumptions tap <C T m <C r s are fulfilled with high probability. 
Assume the neuron is set with a "WSS initial condition", so that {s m , Y m } are also WSS 
processes (or the initial condition time is given at m — — oo instead of m = 0). For other 
initial conditions, we get a transient response, similar to that described in [35) . Since in 
this case the processes {s m } , {Y m } are WSS, they have constant means (s m ) = s* and 
(Y m ) — P* an d possess power spectral densities. 



Linearization near a fixed point In order to study the long term dynamics, we 



linearize the dynamic equations. To do so denote T m = T m — TL Y m = Y m — p*, s m = 



s m — s*,w = Vpap| s . and explicitly solve the system in Eqs. 20 21 , given the following 
assumption. 

Assumption 2: s m is small enough, so that with high probability |s m | <C |s*| (component- 
wise) and |w T s m | > |s^ ( VVpap| s ») s m | . 

Given assumption 1, this essentially means that s* is a stable fixed point of the system 



(Eqs. 20 21), and stochastic fluctuations around it are small, compared to the size of 
the region {s|pap ( s m) 7^ 0, 1} (usually determined by the noise level of the rapid system 
{V, r}, Information SI, section D). Also note that it is guaranteed that, for the "criti- 
cal" stimulation parameters (for which is near 0, 1), the fixed point is stable, at least 
marginally (see Information SI, section E). Given assumption 2, we can approximate 

Pap (s m ) w p* + w T s m , (23) 



which allows us to linearize Eq. [2TJ This essentially means that the components of 
s m determine the neuronal response linearly, with the components of w serving as the 
effective weights (related to the relevant conductances in the original CBMs). Next, we 



wish to linearize Eq. 20 In order to simplify the resulting expressions, we make one last 
assumption. 

Assumption 3: The average kinetic rates are insensitive to the value of s m , namely, for 
all components we require |da±.o (s m ) /ds m \ a± y o (s m ) . 

Consequently, henceforth we approximate the average kinetics rates as constants, as well 
as their constructs, 

A±, = A ±)0 (s*) w A ±i0 (s m ) , 

and similarly for b and D. This approximation results from the typical sigmoid-like 
shape of voltage dependent rate functions and seem to approximately hold in all cases we 
checked (see [35)). However, it is non-essential here, and is made above for mathematical 
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convenience. Using the two latter assumptions, and denoting by I the identity matrix, 
Ai = A + — A_, A* = ((A_ — A ) + p*Ai) tapT" 1 + A , and similar notations for b and 
D we obtain 

s m+ i w T* (A*s* — b*) + (I + T*A*) s m (24) 
+ (A (s* + s m ) — bo) T m + tap (Ai (s* + s m ) — bi) Y m + n, m 



Taking expectations and using Eqs. 16 and pi we obtain 



— (s m -f-x s m ) ~ T^A^s* T^b*, (25) 

where in the last line we neglected a second order term TapW T (s m s^) Ai. Therefore 
s* = A^b* and we can find implicitly from = pap (s* (p*)) for a given T*. We write 
the explicit solution of this equation as 

P. = </CT,). (26) 



Next, using |s m | <C |s*|, Eq. 25 and defining F = I + T*A*,d = AqS* — bo and a = 



tap (AjS* — bi) we can, approximate Eq. 24 



as 



s m+ i = Fs m + df m + aY m + n m , (27) 

which, together with 

Y m = w T s m + e m , (28) 

yields a simple linear state space representation with T m as the input, s m as the state, Y m 
as the output and two uncorrelated white noise sources with variances of 

E n ^(n m O = TJX, (29) 

° 2 e = (e 2 m ) « P*-Pl, (30) 

where in the last line we neglected the second order term w T (s m s^ w. It is straightfor- 
ward to verify also that (T m n m ) = 0, {T m e m ) = 0. 



Input-output statistics 

Eq. [26] provides a deterministic memoryless non-linear relation between T* and - the 



mean statistics of the input and the output. Moreover, the system described in Eqs. 27 



30 



gives a dynamic, linear and stochastic relation between T m and Y m . This allows the 
calculation of all second order statistics in the system, e.g., the auto-covariance i? xy (k) = 
(x m y^ +fe ), or, equivalently, the power spectral densities S xy (e? u ) = J2T=-oo -^W (^) e ~^ k 
(also S x = S xx ). As is standard in the theory of discrete time linear systems |6l| , 
a particularly simple description emerges using the z-transform x (z) = Yl'kL-oo x kZ* • 
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First, we denote the transfer functions H D (z) = (Iz — F) 1 (Open loop) and H c (z) = 
(iz — F — aw T ) 1 (Closed loop), and note that by the Sherman Morrison lemma w T H c (z) 

w T H (z) (l — w T H (z) a) 1 . Then, reading off the diagram at bottom of Fig. [IJc, we 
obtain 



s (z) = H c (z) (n (z) + dT (z) + ae (z)J , (31) 

Y(z) = w T H c (^) (n(z) + df (z)+ae(z)) +e(z) . (32) 
Using standard techniques ( (61]) we find 

S s (e^) = H c (e^) (S n + aa T a 2 + dd T S T (e*")) Hj (e~^) , (33) 

Syt (e JW ) = w T H c dS T (e*") , (34) 

S Y (e^) = w T H c (e^) (S n + dd T 5 T (e^))Hj(e-^)w (35) 
+ ( T c 2 |l + w T H c (e iw )a| 2 . 

For low frequencies it is sometimes more convenient to use the continuous-time versions 
of the PSDs, (/) = T^S^y (e J,w ) aJ=27r ^ Tt for / <C T" 1 , which are approximated by 

S s (f) = (A, + aw T T- 1 -27r/j) _1 (D*+T- 1 aa T a 2 + dd T S r (/)) (36) 

• (A, + aw T r- 1 + 2vr/j)" 1 , 

Syt (/) = w T (A, - 2 7 r/j)- 1 dS T (/) (l - T^w T (2nfj - A,)" 1 a)" 1 , (37) 
S Y (f) = (w T (A*-2nfj)- 1 (B* + dd T S T (f))(Aj + 2nfj)- 1 w + ( T^ (38) 

• |l-T- 1 w T (27r/j-A,)- 1 ar 2 . 



Note that (assuming no zero-pole cancellations) in Eqs. 31 38 all the transfer functions 
and spectra are of order M and have the same poles, which are the A eigenvalues of the 
characteristic polynomial | AI — F — aw T | = in the discrete case, and 

| AI — A* — T~ 1 aw T | = (39) 

in the continuous time case. Since we assumed that the fixed point at s*is stable, necessar- 
ily the solutions for these equation must fulfill |A| < 1 in the discrete case and Re [A] < 
in the continuous case, respectively. 

Optimal linear filtering 



The system in Eqs. 27 30 is also useful for obtaining optimal linear estimators of Y m and 



s m . Using spectral factorization |56|, we convert the system to the innovation form 56 

s m+1 = (F + aw T ) s m + df m + Kv m 
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with K = a + FPw (w T Pw + cr?) , where P is the solution of 

P = FPF T - (w T Pw + al)' 1 FPww T PF T + S n , 

derived from the general discrete-time algebraic Riccati equation, where scalar 
white noise variable with zero mean and variance a% , with cr^ = w T Pw + a\. This 
innovation form is directly related to the optimal linear filter - the Kalman filter. This 
is the case since K is the Kalman gain, while a 2 v and P are equal to the Mean Square 
Errors (MSEs) of a Kalman filter estimation of Y m and s m , respectively, in the case of 
zero input j56j. Moreover, the z transform of this system 

Y 0) = F signa l 0) f (z) + H noise {Z) V (z) , (40) 

with ifsignai {z) = w T H c (z) d and H noisc (z) = w T H c (z) K + 1, describes more concisely 
the contributions of the input and noise to the output as an ARMAx model [62] . 



Examples and simulations 

Having presented the general formalism in previous sections we consider next several 
concrete examples, and show that they all fall into the general formalism developed. We 
begin with the 'HHS' model ( |35|, and see Information S2 for exact parameters) which 
augments the classic HH model with an additional slow inactivation process of the sodium 



conductance 58,63 . The HHS model includes the uncoupled stochastic Hodgkin-Huxley 



(HH) model equations |3||64] 

CV = g Na sm 3 h(E Na -V)+g K n 4 (E K -V)+g L (E L -V) + I(t) (41) 
r = a r (V)(l—r) — /3 r (V)r + a r (V,r)€ r , for r = m,n, h (42) 

with the additional kinetic equation for slow sodium inactivation 



s = 5 (V) (1 - s) - 7 (V) s + a s (V, r) £ s , (43) 

where V is the membrane voltage, / (t) is the input current, m, n and h are ion channel 
"gating variables", a r (V) , /3 r (V) , 5 (V) , and 7 ( V) are the voltage dependent kinetic rates 
of these gating variables, C is the membrane's capacitance, Ek, E^ a and El are ionic 
reversal potentials, and gK,9Na and §l are ionic conductances. Furthermore, a r (V,r) = 
y/N-i(a r (V)(l-r) + l3 r (V)r) and a s (V, s) = ^/N- 1 (5 (V) (1 - s) + 7 (V) s). 



Map parameters Comparing with Eq. 13 we find A(V) = — 7 (V) — 5(V), b(V) = 
S(V), B(V,s) = cr s (V,s). Following the reduction technique described in section 'The 
excitability map', we numerically find the average rates 7-^0 (s) and 5±$ (s) (as in |35| , 
where we used an H/M/L notation instead of the current notation +/ — /0 here), Tap 
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and pap (s) (see Information SI, section D). From Eq. 26 we obtain (Fig [2]A.), w and 
s* for a given T*. The average inactivation rate at steady state is denoted by 



7* - {v*l+ (s*) - (1 - p*) 7- (s*)) r A pT„ 1 + (l - r AP T^ x ) 7 (s*) , 

and similarly for the recovery rate 5*. And so, A* = —7* — 5*, b t . = —5*, s* = £*/ (7* + 5*), 
d = (7A - 7(A) / (7* + £*) and = N' 1 (5*7*/ (7* + 5*)). Denoting = 7+ (s*) - 
7_ (s*) and similarly for 81, we can also obtain a = tap (7*^1 — 7i<5*) / (7* + 5*)- Note 
that s*, A*, 6*, d, and a are scalar now, and so are not boldfaced, as was before. 



Power spectral densities Substituting in Eq. 35 for example, assuming periodic input 
T m gives 



Sv(f) 



w 



(D, + £S T (f)) + T^ e ((2nf) 2 + A*) 
(2tt/) 2 + (A*+%- l wa) 2 



(44) 



Note that this is the shape of a high pass filter when St (/) = (in which case the input 
to the system is purely its internal white noise), as can be seen in Fig [2 3 (top), where we 
compare this analytic result with the simulation and the full CBM (41 43 ) and a simulation 
of the excitability map. Since it follows in this case from (40 ) that Sy (/) = \H no i SC (/) | 2 a%, 
this indicates that H no i se (f) is also a high pass filter. In contrast, S s (f) has the shape of 
low pass filter (Fig|2]B , bottom). Similarly, 



Svrif) = — ™ d . , a M/), ( 4 5) 

2tt fj - (A, + T- l wa) 

has the shape of a low pass filter. Since Syrif) — H signa i (f) St (/), this indicates 
that -f^signai (/) is a low pass filter. We compare this result with a simulation of the full 
CBM (Fig|2p). Since the estimation of the cross-spectrum is known to be noisy ( [65] 
p. 321, and see Fig. S2), we only show a part with a relatively low estimation noise level. 



Optimal linear filtering Using the Kalman filter for the HHS model, as derived in sec- 
tion 'Optimal linear filtering', we can estimate the internal state s m given {{Tk}™~Q , {^j^lTo 1 } 
quite well, with low MSE (Fig. |2p). Also, we tested the Kalman predictor estimate for 
the spiking event Y m given {{Tjt}^, , {^fc}fc=o }, at various noise levels. We compared 
the predictions with the trivial "mean" estimate Y^ 1 = p* (a lower bound on perfor- 
mance), an "oracle" estimate Y^ 1 = pap (%) (an upper bound on performance), and 
two "black box" models (available in the Matlab system identification toolbox): AR- 
MAx(l,l,l) and a 1st order State Space model (see [62] - other black box models tested 
performed similarly). We measured the predictors' performance using the error probabil- 
ity Pcrror — P {Y m ^ I (Y^ > 1/2)), instead of the MSE, since Y m is a binary variable; 
here X() is the indicator function. As can be seen in Fig. [2jZ> , for N = 10 4 and N = 10 6 , 
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all estimators Y^ st perform equally (except for state space model) - poorly yet optimally, 
since -P e ™ e a r n = P°^ c , which indicates that the AP generation noise level is high. Note 
however, that s can be accurately estimated even in this case (Fig. |2p). For iV > 10 s , the 
noise level decreases, and all estimators show improvement over the mean predictor as iV 
increases. The black box model performances are always close to the oracle bound (and 
so, perform optimally), while the model based Kalman filter is slightly above this bound 
at iV = 10 10 . This is not surprising since for iV = 10 10 the fixed point s* looses stability 
(due to a flip bifurcation) and a limit cycle becomes the new stable attractor (as in the 
deterministic model [35]), which invalidates assumption 2, and so our linear model can 
become inaccurate. At N = 10 12 the dynamics become almost completely deterministic, 
and so the estimation error goes to zero for all non-trivial predictors - even our possibly 
inaccurate model. 

Since our model for s m , Y m is linear, it is easy for linear black box models to perform 
accurate system identification. Therefore, it is not surprising that the black box estimators 
manage to perform well as the model-based Kalman filter. In fact, since that the black-box 
models remain close to the oracle bound even when our linear model might be inaccurate 
hints that a different linear model might be derived there. However, the main issue to 
note is that the black box models cannot estimate accurately the state space structure 
(F, d, w, ...) or s m , unless some information is given on the state-space structure, in which 



case they become grey-box models. In other words, our dynamic model, given by Eqs. 24 



and 26, allows us to predict the internal state variables and AP sequences as accurately as 
any predictor, while maintaining a compact and transparent biophysical interpretation. 
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Figure 2. Comparisons of simulations with the mathematical results, for the stochastic 



HHS model. A Eq. 26 compared with a simulation of the full CBM (Eqs. [IT 13), for 
different currents (/ st im = 7.5, 7.7, 7.9, 8.1, 8.3 /j,A from bottom to top). B The spectra 
Sy (/) and S s (f) resulting from periodic stimulation, compared with a simulation of the 



full CBM and the excitability map (Eqs. 20 21). C Kalman filter estimation of s, for 
the same CBM simulation as in B. D Amplitude and phase of the cross-spectrum 
Syr if) for Poisson stimulation. Note that the frequency range was cut due to spectral 
estimation noise (see Fig. S2) E The probability of error in the estimation of Y m , for 
various values of N, and various estimators, as described in section 'Examples and 
simulations'. 



Other models Next we examine additional CBMs (see model details in Information 
S2) and demonstrate that our analytical expressions hold in the following scenarios: (1) 
when the kinetics of the neuron are extended to arbitrarily slow timescales, (2) when 
assumptions 1 and 2 break down, (3) when the rapid and slow kinetics are coupled (in 
contrast to Eqs. 12 13), (4) when "physiological" synaptic inputs are used. 



19 



First, we test whether or not the model can be extended to arbitrarily slow timescales. 
We added to the HHS model four types of slow sodium inactivation processes with in- 
creasingly slower kinetics and smaller channel numbers. In the first case, those processes 
were added additively (as different currents), so s was replaced with ^\ s; in the voltage 
equation |4l| This model was denoted 'HHMS' (HH with Many Sodium slow inactivation 
processes). In the second case, those processes were added in a multiplicative manner (as 
different processes affecting the same channel, in the uncoupled approximation), so s was 



replaced with Yii s i m the voltage equation 41 We denote this model as 'Multiplicative 



HHMS'. In both case, our analytical approximations seemed to hold quite well. For ex- 
ample, the approximated Sy (f) (Eq. [38]) corresponded quite well with the full numerical 
simulation of the stochastic CBM (Fig. |4t3 and D, respectively). 

Next, to test the limits of our assumptions we extended the HHS model (HHSIP 
model, 1 35 1) and added a potassium inactivation current which had faster kinetics (so 
r s Rj 5 Hz). So if T" 1 = 10 Hz, we get T* rs 0.5r s , so the timescale separation assumption 
1 is not strictly valid here. Also, for certain parameter values we get a limit cycle in the 
dynamics of s m , so the fixed point assumption 2 becomes incorrect. However, it seems 
that our approximations still follow the full numerical simulation of the stochastic CBM: 
for at various stimulation frequencies T" 1 and currents J stim (Fig. [3R), for Sy (f) at 
T" 1 = 10 Hz when assumption 1 breaks (Fig. [3^3, top), for Sy (f) at T" 1 = 30 Hz when 
assumption 2 breaks (near a Hopf bifurcation) and a limit cycle begins to form (see Fig. 
[3^3, bottom), and for state estimation of S\ using a Kalman filter, again at T" 1 = 10 Hz 
(Fig. |3p). The only discrepancy seemed to appear in the limit cycle case, where the 
frequency of the limit cycle "sharpens" the peak in Sy (f) (Fig. |3]B, bottom). 

Next, to challenge the approximation even more, we added to the HHSIP model four 
types of sodium currents with increasingly slower kinetics and fewer channels, similarly 
to the HHMS model (so this is the 'HHMSIP' model). This significantly increased the 
variance of the dynamic noise n m , rendering the dynamics more "noisy". These random 
fluctuations in s m (Fig. f3fe) are of similar magnitude to the width of the threshold (non- 
saturated) region in pap ( s m) (similar to Fig. SI). This renders the fixed point assumption 
(2) inaccurate, since now the linear approximation pap ( s m) ~ P* + w T s m breaks down 
most of the time. However, even in this case, the approximations seem to hold quite well 
with simulations of the full stochastic CBM model (Fig. [3JD-F). 

In Fig. [4] A we used a coupled version of the HHS model ('coupled HHS' model), 
in which the equations for r and s in the stochastic CBM are tangled together, and not 
separated as we assumed in Eqs. 12 13 Even in this case, our approximations seemed to 
hold well. 

Finally, in Fig. |4p, we extend the HHS model so that the stimulations are not given 
directly, but through a synapse. We used the biophysical Tsodyks-Markram model |66| of 
a synapse with short-term depression, with added stochasticity ('HHSTM' model). This 
also seemed to work well. 
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Figure 3. Comparing mathematical results with CBM simulation when the 
assumptions fail to hold. In the HHSIP model (HHS with potassium inactivation) 
we plot A (T^ 1 ) for different currents (J s ti m = 7.5, 7.7, 7.9, 8.1, 8.3 /j,A from bottom to 
top). B Sy {f) for two values of T*. Upper figure shows the case when T* w 0.5r s so the 
timescale separation assumption breaks down. In the lower figure the parameters are 
close to a Hopf bifurcation where a limit cycle is formed so the fixed point assumption 
breaks down, so the estimation of limit cycle frequency component is less accurate. C 
The estimation of s\ for T^ 1 rs 30 Hz is even better than in the HHS case. Similarly we 
plot the results of the HHMSIP model (HHSIP with many additional slow sodium 
inactivation kinetics) in (D-F), which has considerable more noise in the slow kinetics, 
and so even larger fluctuations (which further invalidates the fixed point assumption). 
See Information S2 for various model details. 
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Figure 4. Comparing mathematical results (green) with CBM simulation 
(blue) for various models. A Coupled HHS (HHS coupled slow and rapid kinetics) 
B HHMS (HHS with many additional slow sodium inactivation kinetics) C HHSTM 
(HHS with a synapse) D Multiplicative HHMS (variant of HHMS). See section 
Information S2 for various model details. 

Discussion 

Any model is "a simplification and an idealization, and consequently a falsification" (67). 
What details can be 'idealized away' when constructing and using biophysical neuron 
models? To formalize this issue we must properly define: (1) What constitutes the class 
of "excitable neuron" models. (2) The relevant questions that models from that class 
should answer. 

Historically, the answer to the first question was essentially "the set of all models gen- 
erated by electrophysiologists to model neuronal behavior". Such an answer was perhaps 
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satisfying when only a few models of excitable membranes (e.g., the HH model |3|) ex- 
isted. However, since that time the class of allowable models has significantly expanded 



and is still growing rapidly with the discovery of new biophysical processes 27-29]. There- 
fore, if any generic results are to be derived regarding this set, we must first find a more 
useful definition. In this work we suggest a general definition for what constitutes an 
"excitable element" (section 'Stochastic conductance based models'), which covers most 
of the models for excitable point neurons usually used in the literature. 

The answer to the second question depends on the level of investigation. Commonly, 
neurons are explored experimentally by injecting current and measuring the membrane 
voltage. Therefore, relating to experiments, modelers have often tried to use the CBMs to 
find an Input-Output (I/O) relation from current to voltage. This approach is generally 
useful since it is "close to the measured data" and can be used to fit biophysical parameters 
(68], but has two main disadvantages. First, in most cases, due to the complexity of 
CBMs, it can only be studied using a numerical simulation of the model. However, many 
parameters are often unknown, and the simulation results can change unexpectedly due 
to changes in those unknown parameters. Due to the large volume of the parameter 
space and its complexity (e.g., different I/O relations are given by non-convex parameter 
sets (53]), it is hard to perform a proper numerical scan of that space. Therefore, in many 
cases the relevance of such simulation results is not clear. Second, according to the main 
dogma of neuroscience, neurons communicate using Action Potentials (AP), a.k.a. "spikes". 
Therefore, if we are to use a neuronal model to understand how different neurons interact 
in a physiological setting, a current to voltage I/O might be less useful, since it is less 
relevant functionally. And so, to answer the second question, it seems reasonable to require 
that in a "good neuronal model" we should be able to accurately characterize a "spiking 
I/O" relating input spikes to output spikes, and to understand how this I/O changes when 
biophysical parameters are varied. In this work we have found precise conditions on the 
stimulus and model parameters under which CBM models can be transformed into such 
"good models", and show how such a transformation can be executed. 

Analytical results More concretely, we established a formalism through which any 
"excitable element" under a sparse stimulation protocol can be analyzed (section 'General 
formalism'). Applying this formalism to a general stochastic CBM (section 'Stochas- 
tic conductance based models') with arbitrarily complex dynamics, we transformed the 
continuous time mathematical description into a much simpler "excitability map" (Eqs. 



20 21), which has far fewer parameters, is easily interpretable in functional terms, and is 
much more efficient in numerical simulations. This transformation is guaranteed to be 
accurate given a timescale separation assumption (1). For a stationary input, we used a 



linear approximation for the steady state map dynamics (Eqs. 27 28), which is guaran 



teed to be accurate given assumptions 2 and 3. This allows us to write expressions for the 



mean firing rate (Eq. 26) and all second order statistics (Eqs. 33 38), and to construct 



optimal filters for the neuronal response and state (section 'Optimal linear filtering' 



23 



Numerical results We tested the mathematical expressions numerically (Figs. |2j[3]and 
[2J for various parameters and models, and found a good match, even in cases which did 
not strictly obey our assumptions. For example, we show that our results hold reasonably 
well even if assumption 1 and 2 do not hold (Fig|3|. Additionally, although we did not 
explicitly consider synaptical input, it is straightforward to extend the model to do so: 
in Fig. |4p, we added to the HHS model a synapse undergoing short term plasticity 
through synaptic depression (based on [66]), and obtained similar results. Therefore, 
it should be straightforward to include synaptic dynamics (including STDP) and also a 
passive dendritic tree in the framework. However, including more than a single excitable 
compartment might be somewhat more involved, since the neuron no longer possesses an 
"all-or-none" response as assumed here, since some of the compartments can spike while 
others do not (see also |35|). 



Further Generalizations Additional straightforward generalizations are outside the 
scope of this work. First, the use of heterogeneous ("weighted") input spikes (here we 
assumed for simplicity they were all identical). Second, the relaxation of the temporal 
sparsity assumption on the input, tap <C T*. Lastly, the inclusion of ion-channel interac- 
tions at the population level. The existence of such effects in CBMs is currently debatable 
at the fast timescales ( |69| and Brief Communications arising), but are known to exist 
through different cellular mechanisms that might affect excitability (e.g., gene regulation 
networks). Such effects might render channel dynamics non-linear |60|[70] . This might 
also happen if assumptions 1 or 2 become grossly invalid. In both cases, the response of 
the neuron can be non-linear, and reflect a more "digital" kind of computation (e.g., |71|). 
However, in some cases, it might still be possible that a linear model can capture well the 
general model dynamics (Eq. YA Eq. 8 ) near some steady state attractor, if T <C T* (note 
that this was not originally required for linearity), or we might need to use something 
similar to a phase response curve 72 73). 



Estimation The hypothesis relating to the linearity of the residual stochastic process 
is corroborated by the good performance of linear black-box predictors (Fig. [2|3) when 
assumption 2 breaks down. In any case, the relevance of the model suggested here to a 
realistic neuron, can be directly tested experimentally through the linearity of the neural 
response (e.g,. by using sinusoidal T m ). If the neuron responds linearly, then linear black- 
box predictors (such as ARMAx-based predictors) should perform optimally (no need for 



the more general Box- Jenkins-based predictor |62], as indicated by Eq. 40). However, as 
indicated by Fig. [2^3, in some cases (e.g., N = 10 6 ), optimal may not be very different 
from useless due to the high level of stochasticity. Also, even in such cases the neuronal 
state can be estimated well (Fig. |2jC). Both facts indicate that the intrinsic AP generation 
mechanism in the neuron is rather noisy, so simulation of predictors based on deterministic 
CBMs (e.g., (24)) might be misleading. 
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Relation to previous work To our knowledge, no previous work examined the re- 
sponse of general stochastic CBMs to temporally sparse ("spike-like") input. However, 
in 35 , this was done for general deterministic CBMs. Besides the inclusion of stochastic 
noise, our model is more general than in (35] also because we allow a general spike train 
stimulation (not just periodic) and for the slow dynamics to be coupled directly (Eq. 13). 
The inclusion of noise in the model actually simplified the analysis, improved significantly 



the accuracy of the mean firing rate approximation (Eq. 26), allowed the calculation of 
the second order statistics as well as the Kalman filter, and improved the similarity of the 
model with the experimental results (as explained in [35]). Naturally, for N — > oo, our 
simulations recover the results in [35] (not shown). We note also that other recent works 
investigated somewhat similar issues in specific (simplified\phenomenological) models. 
For example, spike-like inputs 74,75 , adaptation currents [76], adaptation and noise [25] 
or adaptation and spike-like inputs~[77] . Additionally, map based event-driven simulation 
methods were previously considered in [78], and recently in [79]. However, in the "ex- 
citability map" derived here, the events are input spikes - which considerably simplifies 
the simulation. 



Functionality and Networks The resulting input-output relation we have described 
here is similar to linear analog circuits working at some operating point (s*,p«), deter- 
mined by the input bias (T*). Such a system could be used as an adaptive filter - where 
the mean firing rate tunes the parameters of the linear filter, as well as its spike time 
input-output relation. Note that although the filter is linear when its input is T m and 
its output is Y m , it is not linear when considering the mapping between the input spike 
train to the output spike train. And so, when embedded in a network, the network is not 
necessarily linear. In any case, since different input spikes generally do not "collide" in 
our framework (the tap <C T m assumption) the computations performed by such a filter 
are inherently different from the standard paradigm of neuronal computation j4l], which 
usually involves the summation of different spikes arriving together. Such a summation 
may depend on details of millisecond delays between arriving spikes (a rather standard 
situation [31]), and might be less robust than a stochastic computation based on slow 
kinetics, as presented here. 
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Information SI 

- Mathematical Appendix 



A Noise correlations 

Recall that we defined e m = Y m — pap (x m ). Note that (e m \H m -i, x m) — (e m | x m) = 0, and 
therefore (e m f (9i m -i, x m )) = for any deterministic function /. Also, from its definition, 
e m = / (Hm) f° r some deterministic function /, and so, for n > m 

(^n^m) ( (fin&m \H-n— 1 j x n) ) ( (^n \T~^n— 1 > x ri) 6 m ) 0, 

and similarly for m < n. So e m is a white noise process with (e n e m ) = (e^) 5 mn and 
( e m) = (Pap (xm)) - (Pap ( x m))- We defined n m = Ax m - (Ax m |x m , T m , K m ). Note that 
(n m \'H m ) = (n m |x m , T m , Y m ) = and therefore (n m f (H m )) = for any deterministic 
function /. Note also that from its definition, n m = / ("H m ,x m+ i) for some deterministic 
function /. Since {H m ,x m+ i} C T-L m +k for all k > 1, we obtain for n > m 

( n « n m) = << n nn^|H„)) = <(n n |H n ) n^) = 0, 

and similarly for m < n. So n m is a white noise process, with ^n n n^) = (n m n^ 5 nm . 
Additionally, we note that n m and e m are uncorrelated 

|"H n _i,x n )) ,n>m 
(e n n m ) = < =0. 
I (e n (n m |H m )} ,n<m 



B Derivation of stochastic differential equations 

To facilitate mathematical analysis and efficient numerical simulation, it is sometimes 
desirable to approximate stochastic CBMs using stochastic differential equations (SDE). 
This was initially suggested by [64], but their method suffered from several problems |80| . 
In a recent paper 60 a more general method was derived, which had none of the previous 
problems, and was shown numerically to produce a very accurate approximation of the 
original Markov process description. Here we derive the neuronal SDEs, shown in Eqs. 11- 
13. Using the stochastic CBM notations with Eqs. 11-13 and the method described in |60 
we obtain the following SDE for the dynamics of x^ c ^ 

x (c) = A (c) x (c) + B (c) £ (c) , (46) 
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where the elements of the rate matrix A*^ are defined to be A^- for all i j and 
Ay? = — on the diagonal; £^ is a white noise vector process - meaning it 

has zero mean and auto-covariance (t) (t')) \ = I$c,c'5 {t — if) where I is the 

identity matrix, (•) denotes an ensemble average, S (t) is the Dirac delta function, and 
<5 CjC / = 1 if c = d and otherwise; and B^ c ^ is defined so that each component of £^ 



is associated with a transition pair i ^ j, multiplied by 



/l(c)™(c) _,_ 4^1 



|(c)^(c) 



/iV( c ), 



and appears in the equation for iff and xf with opposite signs. Note that B^ c ^ is not 
necessarily square since it has rows but the number of columns is equal to the number 
of transition pairs. 

Since x^ff denote fractions, Ylk x k = f° r a ^ c - We use this constraint, together 
with the irreducibility of the underlying ion channel process, to reduce by one the di- 



mensionality of 46 (see |81| for details). Defining J to be the I with it last row removed, 
e 4 (0,0,..,l)Vu = (1,1,-, 1) T , G 4 (l-eu T ) J T , A (c) 4 JA^G, B< c > 4 JB^ 
(with x^\ c) replaced by 1 — x\ — Xi--- — Xk( c )-i) an d b^ c ^ = — JA^e (A^ c ^ is invertible), 
we obtain the following equation for smaller state vector y( c ) = Jx^ c ^ (which has only 

- 1 states) 



Ac) 



A (c) y (c) - b + B (c) £ 



(c) 



Furthermore, it can be shown |8l| that A^ c ^ is a strictly stable matrix (all its eigenvalues 
are also eigenvalues of A^ c ^ except its zero eigenvalue, and so have a strictly negative real 
part), and T)^ = B^B^ 7 is positive definite (so all its eigenvalues are real and strictly 
positive). 

Assuming that all channels can be classified as either "rapid" or "slow" (this assumption 
can be relaxed, so that this timescale separation can occur also within a single ion channel 
as we do in the coupled HHS model - see Fig. 4A), we concatenate all vectors related to 

T 



rapid channels r 



(i)' 



(R) 



and to slow channels s 



where 



^(-R+i)' —■>y{R+s) / 

R + S = C . We obtain Eqs. 11-13 by similarly defining b,b',£ and £' together with the 
block matrices 



A' 



/ A W 



V o 







o 





\ 



( A {R+1) 




V o 












\ 



~^(R+S) 



and similarly for B' and B . Note that A' is square with M' = J2c =1 - R rows while 
A is square with M = J2f=R+i ~ & r °ws. 
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C Details of numerical simulations 



MATLAB code available on the ModelDB website, with accession number 144993. In all 
the numerical simulations of the full stochastic CBM we used the SDE Eqs. 11-13. We 
used first order Euler-Maruyama integration with a time step of dt = 5 fisec (quantitative 
results were verified also at dt = 0.5 /isec). Each stimulation pulse was given as a square 
pulse with a width of i stim = 0.5 ms and amplitude / s tim- The results are not affected 
qualitatively by our choice of a square pulse shape. We define an AP to have occurred if, 
after the stimulation pulse was given, the measured voltage has crossed some threshold Vth 
(we use Vth = — 10 mV in all cases). In all cases where direct stimulation is given, unless 
stated otherwise , we use periodic stimulation with J st i m = 7.9/iA and T* = 50 msec. Note 
that due to parameters, no APs are spontaneously generated, as seen experimentally |3l) . 

D Calculation of £>ap ( s ) 

We numerically calculated pap (s) from a stochastic CBM by setting e = and disabling 
all the slow kinetics in the model (so we only use Eqs. 11-12). Then, for every value of s 
we simulated this "semi-frozen" model numerically by first allowing r to relax to a steady 
state and then giving a stimulation pulse with amplitude istim- We repeat this procedures 
200 times for each s, and calculate pap (s) as the fraction of simulations that produced 
an AP. A few comments are in order: (1) In some cases (e.g., the HHMS model) we can 
use a shortcut and calculate pap (s) based on previous results. For example, suppose we 
know the probability function Pap (s) for some model with a scalar s and we make the 
substitution s — h (s) where the components of s represent independent and uncoupled 
channel types |60] - then pap (s) = pAp(h(s)) in the new model. (2) The timescale 
separation assumption tap <C T m <C r s implies that all the properties of the generated 
AP (amplitude, latency etc.) maintain similar causality relations with s m as does Y m , so 
we can find their distribution using the same simulation we used to find pap ( s )j similarly 
to the approach taken to compute L (s)in the deterministic setting |35|. (3) Numerical 
results (Fig. SI) suggest that we can generally write 

PAP(S) = $(£(8)/^) (47) 

where $ is the cumulative distribution function of the normal distribution, E (s) is 
some "excitability function" (as defined in (35), so pap (s) = 0.5 on the threshold O = 

1 /2 

{s\E(s) = 0}), and N' 1 ', the ' 'noisiness" of the rapid sub-system, directly affects the 
slope of pap (s) (Fig. SID, bottom). Also, as explained in |35|, E (s) is usually monotonic 
in each component separately and increasing in J st i m (Fig. SIC, top) - which could be 
considered as just another component of s which has zero rates. Note that the J st i m de- 
pendency of pap ( s ) is commonly measured experimentally (averaged on the values of s 
in the experiment), so this can be used to get a crude lower bound on A^ r . 
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Eq. 47 can be explained heuristically. Recall that the fast variables {V (t) , r (t)} follow 
stereotypically either the "AP response" (Y m = 1) or the "no-AP response" (Y m = 0), 
then equilibrate rapidly to some quasi-stationary distribution q(V, r|s), which becomes 
q (V, r) given assumption 3. Since the voltage response is stereotypical given Y m , the rapid 
ion channels are effectively independent during the response. Also, unless the channel 
numbers are very low (N < 10 2 ), the white noise £ in Eq. 12 can be safely approximated 



as Gaussian 60 , 80 1 . Finally, at the rest voltage (V mst ), both V and r display only small 
fluctuations. Combining all the above facts, we can approximate r = r* + r and V = 
V Test +cJi where f is a "small" zero mean Gaussian variable and C\ is some constant vector. 
Given a sparse spike stimulation, an AP will be generated if and only if E (V, r, s) + X > 
where E (V, r, s) is some excitability function and X is some random variable (possibly 
dependent on V, r and s). The AP condition then becomes 

< E (V rcst + c T f , r* + r, s) + X 
« £(V rest ,r*,s) + cJr + X 



If Var (cjr) 3> Var (X) then we get 47 since cjr is a Gaussian variable as a sum of 
Gaussian variables, with Var (cjr) oc Nr^^ as in (Fig. SID, bottom). 



E Initial condition far from steady state 

Transients What happens when the input intervals T m constitute a WSS process (as 
in section 'Input-output relations for a stationary input'), but s m is initialized at a value 
which is far from the steady state? Suppose that this is the case, and that in addition, all 
our assumptions hold (except for assumption 2). As demonstrated numerically (Fig SI), 
and hinted experimentally ( (82] Fig. 1A&B), it is common for pj^p (s) to be practically 
equal to or 1 for most values of s - except in a small region. Therefore, if we denote 
E + = {s\p (s) = 1} , E- = {s\p (s) = 0} , E = {s\p (s) 7^ 0, 1}, then as long as s m is not 
in Eq, we obtain 

, v = f ((A+ - A ) r r + % A ) «s m ) - s+ ) , if (s m > G E + 

{m+1 m) \((A_-A )r r + T,Ao)((s m )-s M ) , if (s m ) G E_ 

with (s m ) — > = ((A ± — A ) r r + T*A ) _1 ((b ± — b ) r r + T^bo) in each case (since 
(A± — Aq) r r + T^Ao are contracting matrices). 

Steady state Eventually s m reaches some stable steady state behavior. There are 
several different possible steady state modes, depending on sj^ (which are affected by 
Jstim and T*). If both s^ ^ Eq then steady state behavior can be found by the following 
self-consistency arguments (similar to those used in (35)). 
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1. Responsive: If s+ G E + , s ro G E + then s m stabilizes at s+ , so Y m = 1 for all m. 

2. Unresponsive: If s+ G E-, s ro G EL then s m stabilizes at s^, so Y m = for all m. 

3. Intermittent: If s ro G E + , s+ G EL then s m stabilizes on a fixed point in E , or on 
an attractor (e.g limit cycle) intersecting E . 

4. Bi-stable: If s+ G E + , G E_ then, depending on the initial conditions, s m might 
stabilize either on s+ (and then Y m = 1 for all m), on (and then Y m = for all 
m), or (in more pathological cases) on some attractor that intersects Eq. 

And so, the steady state neuronal response {Y m } can be a non-constant series only in the 
intermittent mode, if s+ G E , if G E or in some bi-stable cases (note that the bi- 
stable case was not observed in cortical neurons (35j). In these cases, {Y m } becomes a true 
"stochastic" process, which is coupled with {s m }. Note also, that in both the responsive 
and unresponsive modes s m remains near a stable fixed point. Therefore, if we change 
stimulation parameters continuously (e.g., T* or J s ti m ) and reach the intermittent mode 
at some critical parameter, then the fixed point is still stable (or marginally stable). 
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F Additional Figures 




Figure SI. Fitting of p^ P (s) = $ ((s — a) /&) in the stochastic HHS model. A 

Fitting of pap (s) f° r various values of J st i m . B Fitting shows that a is linearly decreasing 
in Jstim- C Fitting of pap (s) for various values of N. D Fitting shows that b oc 1/y/N. 
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Figure S2. Estimation noise in the cross-power spectral density. When 
estimating cross-spectra, estimation noise level can be quite high (ocl /coherence, 
according to j65], p. 321). To estimate the level of this noise in Fig. [2p, we added 

S Y f (f) where| T m j is a shuffled version of {T m }. Only when the estimated Syr (/) is 

above S Y f (f), is its estimation valid. Therefore, in figure [2jZ> we show only this region 
(left of dashed black line), where estimation is valid. 
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Supporting Information S2 - 
Model Details and Parameters 



As argued in the introductory section, a main asset of the present approach is its ap- 
plicability to a broad range of models of various degrees of complexity and realism. In 
this section we describe the models used in the paper and give their parameters. These 
models have either been studied in the literature or are extensions of such models, which 
are meant to explore the limit for the validity of our analytic approximations. In all cases 
where direct stimulation is given, unless stated otherwise, we use periodic stimulation 
with I stim = 7.9/xA and T* = 50 msec. 



The HHS model 



vation 58, 63 



The HHS model combines the Ho dgkin- Huxley equations |3j with slow sodium inacti- 
The model equations (35], which employ the uncoupled stochastic noise 



approximation, are 



CV = g Na sm 3 h(E Na -V)+g K n i (E K -V)+g L (E L -V) + I(t) 
m = [a m (V) (1 - m) - m (V) m] + \J N h x (a h (V) (V) (1 - h) + J h (V) h)U 

n = [a n (V) (1 - n) - (3 n (V) n] + \J N^ 1 (a h (V) (V) (1 - h) + J h (V) h)t n 

h = [a h (V) (1-h)- p h (V) h] + \J N^ 1 (a h (V) (V) (1 - h) + (3 h (V) h)£ h 

s = 6 (V) {l-s)- 1 (V)s+ VA-i (6 (V) (1 - s) + 7 (V) s)£ s 

Most of the parameters are given their original values (as in |3|[63]): 

V Na = 50 mV, V K = -77 mV, V L = -54 mV, 
g Na = 120 (kQ ■ cm 2 )' 1 , g K = 36 (jfefi ■ cm 2 )' 1 , g L = 0.3 (kQ ■ cm 2 )' 1 , 

a n (V) = i ^ilff-fL) kHz, MV) = 0.125 ■ c^ym kHZ; 

a m (V) = kHz > PJY) = 4 • e-^+ 65 )/ 18 kHz, 

a h (V) = 0.07 • e -(^+65)/20 kHZ; ^(y) = ( e -o.i.(v+35) + kHZ; 

where in all the rate functions V is used in units of mV. In order to obtain the specific 
spike shape and the latency transients observed in cortical neurons, some of the parameters 
were modified to 

C m = 0.5 /iF/cm 2 , = 2 
7 (V)=0.51-(e- a3 -( y+17 ) + l)- 1 Hz , 5(V) = 0.05e^ v+8 ^ 30 Hz 
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Estimates of channel number greatly vary |35|. For simplicity, we chose N = N n = Nj 



N m = N s , and unless stated otherwise, we chose, by default N = 10 , as in 35 



The Coupled HHS model 

The coupled version of the HHS model has the same parameters as the uncoupled version, 
and a similar voltage equation 



CV = g Na s m h {E Na -V) + g K n {E K -V)+g L (E L — V) + I (t) 

where the variables the n Q and So^-o^o describe the respective fraction of potassium and 
sodium channels residing in the "open" state. To obtain the coupled model equations, 
we need to assume something about the structure of the ion channels. The original 
assumption by Hodgkin and Huxley was that the channel subunits (e.g., m,n and h) are 
independent. Over the years, it became apparent that this assumption is inaccurate, and 
the sodium channel kinetics subunits are, in fact, not independent |83|. However, it is 
not yet clear how the slow sodium inactivation is coupled to the rapid channel kinetics 
(e.g. 



84 85]), so we nevertheless used the original naive HH model assumption that the 
subunits are independent. In that case the potassium channel structure is given by (for 
brevity, the voltage dependence of the rates is henceforth ignored for this model) 

Aa n 3a n 2a n a n 

n o ^ n i ^ n 2 ^ n 3 ^ n i 
f3 n 2f3 n 3f3 n Af3 n 

while for the sodium channel is described by 
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3a m 
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s Q m 3 hi 
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3a m 


2a m 
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Pm 


2 Pm 


3/3 m 








3a m 


2o m 






Sim hi 




sxTUihi ^ S\m2hx 




sim 3 hi 




Pm 




3(3 m 





Here if transition rates are indicated between two boxed regions, it is to be understood 
that the same rates are used between all corresponding states in boxed regions. The 
corresponding 32 SDEs are derived using the method described in (60] (or 30 equations if 
we use the reduction technique described in ). In this model we used J st i m = 8.3/xA. 



The HHSTM model 

In order to investigate the effect of a more "physiological" stimulation, we changed the 
HHS model and added synapses. We used the popular Tsodyks-Markram model for the 
effect of a synapse with short-term-depression on the somatic voltage (the model first 
appeared in |66| and was slightly corrected in (86)). In the model x, y and z are the 
fractions of resources in the recovered, active and inactive states respectively, interacting 
through the system 

y 

/ \ • (49) 

x < — z 



Here the z — > x rate is r~l, the x — > y rate is t^" 1 , and the x — > y rate is UseS (t — t sp ), 
where 5 (•) is the Dirac delta function, and t sp being the pre-synaptic spike arrival time. 
The post-synaptic current is given by I s (t) = AsEy(t). Additionally, we added noise 



to the model using the coupled SDE method (60], assuming that diagram 49 with the 
corresponding rate hint at the underlying Markov kinetic structure, with N = 10 6 . As 
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in 66 Fig. IB, we used Ti n = 800 msec, r rec = 3 msec and Use = 0.67. Additionally, we 
set A$e — 70 fx A to obtain an AP response in our model. 



The HHMS model 

The HHMS model consists of many sodium currents, each with different slow kinetic 
variable. The equations are identical to the HHS model, except g^aS is replaced by 
g Na M~ l YaLi s i> where 

iH = 8 t (V) (1 - Si ) - 7 , (V) Si + ^A i - 1 (5 i (y)(l- Sl )+ 7i (^) i)& 

with ji (V) = 7 (V) e\8i (V) = 8 (V) e* and N (i) = Ne*, where 7 (V) , 8 (V) and N are 
taken from the HHS model. Also, we chose e = 0.2, r] = 1.5, M = 5 

The Multiplicative HHMS model 

The Multiplicative HHMS model is identical to the HHMS model with rj = 1, except that 
(jNaM' 1 YhU Si is replaced with g Na ]\ i=1 Si. 

The HHSIP model 

The HHSIP model equations p5] are identical to the HHS mode equations, except that s 
is renamed to s\ and an Inactivating Potassium current was is to the voltage equation 



Ik = 9Mn 4 s 2 (E K - V) . 

with g M = 0.05gK and 

S2 = 82 (V) (1 - s 2 ) - 7 2 (V) s 2 + y/N-i(6(V)(l-s 2 )+i{V)s 2 ), 
where N S2 = N and 

3 3 e (^+ 35 )/ 15 + e -(V+35)/20 3 3e (y+35)/15 _|_ e -(F+35)/20 
5 * = 1 + e -(V+35)/10 HZ ' ^ ^) = " 1 + e( y + 35)/10 HZ - 

Again, in all the rate functions V is used in mV units. In this model we used J st i m = 8.3/xA 
and T* = 33 msec. 

The HHMSIP model 

The HHMSIP model combines HHSIP and HHMS. Its equations are identical to the 
HHMS model with 77 = 2, except it also contain the Ik current from the HHSIP model. 
In this model we used J st i m = 8.3/xA and T* = 33 msec, unless otherwise specified. 
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Supporting Information S3 - 
Nomenclature 



symbol 


Descrintion 




time 


+ 


SlimUiailOIl llllltjb 


171 


QlSCIclc Lllllc lllQcX COIIcSpOllQlllg, LO i m 


T 

- 1 - in 


T T i 1 n T /~\V nf 1 1 > "l IT 1 Oil /~~\T"~^ 1 Y~\ T ATI rn 1 1 

C m _|_l — C m ^lllttjl-StllilUldtlOIl llrieiVaiJ 


v 

1 m 


OCCUIieilCe OI ail /\r ^erCIiei U OI 1J, drCfc!I Stimulation SpiKt! /7J 


V 


the excitable membrane's voltage 




voltage threshold - an AP has occurred if and only if V crossed V t h upwards 


tap 


upper timescale of an AP 


x(f) 


internal "state" of a general neuron model at time t (V is a component of x) 


x ra 


x {t m ) - the state sampled at time t m 




{{x fc }^ =0 , {T fc }™ =0 , {Y k }™ =0 } ( the history set) 


p A p (x m ) 


P (Y m = l| x m) (the probability to generate an AP given the state) 




the mean of X 


(X|F) 


the conditional mean of X given Y 




Y m ~ Pap ( x m) (AP generation noise) 




x m+i — x m (increments) 




Kronecker's delta 




(Ax m 


x m; T m , Y m 1) 


h (x m , T m ) 


(Ax m 


x m; T m , Y m 0) 



Table 1. General neuron models 
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symbol 


Description 


c 


index of ion channel types 


C 


number of ion channel types 


K (c) 


number of states in c-type ion channel 




number of c-type channels 


A? (V) 


transition rate from state j to state i in a c-type channel 


/(•••) 


f = /(•••) 


7(0 


stimulation current 


^stirn 


amplitude of stimulation pulses 


^stim 


duration of stimulation pulse 


s (f) /r (0 


slow / rapid state variables of neuron (so x = (V, r T , s T ) 


tap 


as defined before. Now also upper timescale of {V, r} 


r s 


lower timescale of s 



Table 2. Stochastic Conductance-based models 
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symbol 


Description 


AW (V) 


Matrix of A% , with A\f = - ^ 


(t) 


white noise vector process for c-type channel 


6(t) 


Dirac's delta function 


T3 (c) ( T T (c) \ 

B l > [y, x. [ ') 


D f a( c ) ( c ) i /l( c ) ( c )\ lAT(r\ 1 i- • \ 

Bij = J [AljXj + Aj/xl j /J\( c > (non-square matrix) 


I 


the identity matrix 


e 


(0,0, ...,1) 


u 


(1,1, 


J 


I with its last row removed 


G 


(I - eu 1 ) J 


A (c) (V) 


JA^G (A( c ) is a strictly stable, invertible, matrix) 


B( c ) (V, x ( c )) 


( \ (c) 

JB^- 1 (with xy (c) replaced by 1 — x\ — x 2 ... — _i) 


D(c) (y, x ( c )) 


B( C )B( C ) T (diffusion Matrix) 


(1/) 


-JA^e 


y (c) (*) 


Jx^ (of size - 1) 


r(t) 


as defined before, now also ((y lyl ^) T , (y^) 7 ) 


8(0 


as defined before, now also ^(y^ +1 ^) T , (y^ +S ^) T ^ 


M' 


Ef-i - R (size of r) 


M 


ESSVi ^ W - S (size of s) 




((ba)) T ,...,(b(-)) T ) 


b(V) 


((b( R+1 )) T ,...,(b^ +5 )) T )' 


£' (*) 


(U {1) Y ,-,U iR) Y\ 

V v / v / J 




(U {R+1) ) ,-,U (R+s) ) 

\\ / \ J J 


A (V) 


Diag A , • • • , A J (diagonal block matrix of arguments) 


A(V) 


Diag(A (K+1) ,--- ,A (H+b) ) 


B' (V, r) 


Diag (BW, • • • ,P) 


B(V,s) 


Diag (b( r+1 ),-- - ,B^+ 5 )) 


D'(y,r) 


B'B" 


D(F,s) 


BB 1 


roo(^) 


(A') -1 b' (steady state at a fixed voltage) 


Soo (V) 


A _1 b (steady state at a fixed voltage) 



Table 3. Derivation of stochastic differential equations (reduced form) 
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symbol 


Description 


T T 


* \J"m) i ^* \J"m) i S 


q{V,r\s) 


stationary distribution or {1/, r} given a fixed s 


a(V) 


some microscopic kinetic rate (A^) 


oc+ (s m ) 


// m+TAP a (V 00) dt (averaged a (V)) during an AP response 


«- (s m ) 


TXp f^ m+TAF a (V (t)) dt (averaged a (V)) during a no-AP response 


«0 (Sm) 


t — T^p ii*"+T A p a ^ ft)) dt (averaged a (V) during a steady state) 


{0, ±} subscript 


replace all rates in expression with the corresponding averaged {0, ±} rates 


Pap (s m ) 


P (Y m 1 s m ) 


N 


min c 



Table 4. The excitability map 



symbol 


Description 


z% 


(z) (also - value at fixed point) 


z 


z — z* (perturbation from fixed point) 


z±,o 


z±,o ( s *) (due to assumption 3) 




z + — z_ 


p* 


(Y m ) (instead of Y # , since this is also mean probability of AP) 


w 


V^ap s * (weights of s near fixed point: p^p (s m ) ~ + w 1 s m ) 


A.* 


((A_ - A ) + p.Ai) tapT" 1 + A 


b* 


((b_ - b ) + p*bi) tapT- 1 + b 


D* 


((D_ - D ) + p»Di) tapT" 1 + D 


g(-) 


P*=g (%) 


F 


I + T*A* 


d 


AoS* - b 


a 


r AP (Ais* - bi) 


S n 


(n m n^ t ) (AP generation noise variance) 


°l 


(e^) (state dynamics noise variance) 



Table 5. Input-output relations for a stationary input 
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symbol 


Description 


R xy (k) 


(x m y^ +fc ) (auto-covariance) 


S xy (e*") 


X]fcL-oo -^xy (k) e~i wk (power spectral density) 


Sxy (/) 


T^S'xy (e JtJ ) w=27r/ - Tji 




5xx 


x (z) 


Xlfel-oo x k z ~ k (two sided z-transform) 


n ( z ) 


(Iz — F)~ (open loop poles) 




(Iz — F — aw T ) 1 (closed loop poles) 



Table 6. Input-output statistics 



symbol 


Description 


P 


FPF T - (w T Pw + a 2 e ) ~ X FPww T PF T + S n 


K 


a + FPw (w T Pw + (jg) 1 (Kalman gain) 




white noise with (v m ) = and variance a 2 


°l 


W 1 PW + <7g 


-^signal \Z) 


w T H c (z) d (input to output transfer function) 


-^noisc (^) 


w 1 H c (z) K + 1 (noise to output transfer function) 



Table 7. Optimal linear filtering 
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